function gen_const_qw(nw, nt, imbalance = 0.0)
    qw = ones(nw, nt)
    for iw = 1:nw
        qw[iw,:] = view(qw, iw, :).+rand(-imbalance:0.1:imbalance)
    end
    return qw
end

# function gen_const_qw(nw, nt)
#     qw = ones(nw, nt)
# end

function gen_rand_qw(nw, nt)
    qw = zeros(nw, nt)
    for i=1:nw
        il = sort(unique(vcat(1,nt,rand(1:nt,7))))
        for k = 1:length(il)-1
            qw[i,il[k]:il[k+1]] .= rand()
        end
    end
    return qw
end

function gen_real_rand_qw(nw, nt;
                          mult = 1.0)
    qw = zeros(nw, nt)
    #Кусочно постоянные режимы
    for i=1:nw
        il = sort(unique(vcat(1,nt,rand(1:nt,Int64(ceil(nt/12))))))
        for k = 1:length(il)-1
            qw[i,il[k]:il[k+1]] .= rand()
        end
    end

    #Отработка нагнеателдьных скважин
    for iw in unique(rand(1:nw, Int64(ceil(nw/2))))
        ia = findfirst(qw[iw].<0.5)
        if !isnothing(ia)
            qw[iw,ia:end] .= -2*qw[iw,ia:end]
        end
    end

    #Отключаем скважины случайным образом
    for iw = 1:nw
        off_i = sort(unique(rand(1:nt,Int64(ceil(nt/12)))))
        ia = findall(diff(off_i).<3)
        deleteat!(off_i,ia.+1)
        qw[iw,off_i].=0
    end
    uf = qw.!=0
    return qw.*mult, uf
end

function gen_periodic(nw, nt; ampl = 1.0, trend = 0.0, level = 1.0, cnt = 6)
    qw = zeros(nw, nt)
    ω = 2*pi/nt*cnt
    for iw = 1:nw
        ϕ = (iw-1)*ω
        qw[iw,:] = make_periodic(1:nt, ampl, trend/nt, level, ω, ϕ)
    end
    return qw
end

make_periodic(t, A, B, C, ω, ϕ) = C .+ A .* sin.(ω .* t .+ ϕ) .+ B .* t

function scaleRateByGeoProp(qw_temp, Vp, kp, he, dt, P0, ds)
    #Согласуем режим работы с геологическими параметрами
    nc = length(kp);
    nw, nt = size(qw_temp)

    сVp = sum(Vp)
    dpv = sum(filter(x->x>0.0, qw_temp[:]))/сVp*dt; #Дренированный пороввый объём
    ipv = sum(filter(x->x>0.0, qw_temp[:]))/сVp*dt; #Прокачанный пороввый объём
    if max(dpv,ipv)/(nt/60)<1
        multQ = (nt/60)/max(dpv,ipv)
    else
        multQ = Inf
    end
    Rw = sqrt((ds*nc/nw)/pi)
    qw_m = mean(filter(x->x!=0.0, qw_temp[:]))
    dP = qw_m/(mean(kp)*mean(he)*10 .*8.64*1e-3)*log(Rw/0.05)
    multP = P0/dP
    mult = round(min(multQ,multP), sigdigits = 2)
    qw_new = qw_temp.*mult
    return qw_new
end

function scaleRateByComp(qw; bnd_comp = 1.0, set_comp = false)
    nw, nt = size(qw)
    comp = get_comp(qw);
    comp = smoothN(comp, 3, 3)
    qw_new = copy(qw)
    for t = 1:nt
        if set_comp
            qw_new[view(qw,:,t).<0.0, t] = qw[view(qw,:,t).<0.0,t]./comp[t].*bnd_comp
        else
            if comp[t] > bnd_comp
                qw_new[view(qw,:,t).<0.0, t] = qw[view(qw,:,t).<0.0,t]./comp[t].*bnd_comp
            end
        end
    end
    return qw_new
end

function get_comp(qw)
    sqwp = sum.(filter.(y->y.>0., eachcol(qw)))
    sqwi =  sum.(filter.(y->y.<0., eachcol(qw)))
    comp = -sqwi./sqwp;
end

function movingaverage(a::Vector, n::Int = 3)
    #Скользящее окно n - кол-во соседей
    len = length(a)
    b = zeros(Base.eltype(a), len)
    delta = ceil(Int64, n/2)
    for i=1:length(a)
        si = max(1, i-delta)
        ei = min(len, i+delta)
        b[i] = sum(view(a,si:ei))/(ei-si+1)
    end
    return b
end

#украденный moving average
# function movingaverage(X::Vector, numofele::Int)
#     BackDelta = div(numofele, 2)
#     ForwardDelta = isodd(numofele) ? div(numofele, 2) : div(numofele, 2) - 1
#     len = length(X)
#     Y = similar(X)
#     for n = 1:len
#         lo = max(1, n - BackDelta)
#         hi = min(len, n + ForwardDelta)
#         Y[n] = mean(X[lo:hi])
#     end
#     return Y
# end

function smoothN(a, it, n=3)
    #сглаживание скользящим окном
    b = similar(a)
    b .= a
    for i in 1:it
        b .= movingaverage(b, n)
    end
    return b
end

function addNoise(val, ampl = 0.05; norm_dist = false)
    nw, nt = size(val)
    rsl = copy(val)
    for iw in 1:nw
        vi = view(val,iw,:)
        if norm_dist
            #rsl[iw,:] .= view(rsl, iw, :) .+ifelse.(vi.!=0, ampl.*randn(nt), 0.0)
            add_noise!(view(rsl, iw, :), ampl)
        else
            noise = range(-ampl,ampl, 11).*mean(abs, vi)
            #rsl[iw,:] .= view(rsl, iw, :) .+ifelse.(vi.!=0, rand(noise,nt), 0.0)
            add_noise!(view(rsl, iw, :), ampl, rand(noise,nt))
        end
    end
    return rsl
end

function add_noise!(x, scl=0.1, rnd = randn(length(x)))
    x .+= mean(x).*ifelse.(x.!=0, scl .* rnd, zero(eltype(x)))
end

function add_noise(x, scl=0.1, rnd = randn(length(x)))
    y = copy(x)
    add_noise!(y, scl, rnd)
    return y
end

function offWellbyinterval!(val, cnt = 1)
    #Выключаем скважины по интервалам, cnt - кол-во повторов отключения
    nw, nt = size(val)
    step = floor(Int64,nt/(nw+1)/cnt)
    for (k,v) in enumerate(zip(Iterators.partition(step+1:nt,step), Iterators.cycle(1:nw)))
        val[v[2],v[1]] .= 0.0
    end
end

function gen_rgm(nw)
    #Создание таблицы всех комбинаций режимов
    nr = 2^nw
    fl = falses(nw,nr)
    for i = 1:nr-1
        tmp = parse.(Bool,split(string(i; base=2),""))
        fl[1:length(tmp),i+1].=tmp[length(tmp):-1:1]
    end
    return fl
end


function wellTypePlase(script, nw = 9)
    if nw == 5
        @info "5 скважин"
        if script == 1
            @info "Классическая пятиточка"
            inj_iw = [3]
            prod_iw = [1,2,4,5]
        elseif script == 2
            @info "Диагональ добывающие"
            inj_iw = [2,4]
            prod_iw = [1,3,5]
        elseif script == 3
            @info "Диагональ нагнетательные"
            inj_iw = [1,3,5]
            prod_iw = [2, 4]
        else
            @info "Обращённая классическая пятиточка"
            inj_iw = [1,2,4,5]
            prod_iw = [3]
        end
    elseif nw == 9
        @info "9 скважин"
        if script == 1
            @info "Девятиточка, 3 нагнетательные V обр."
            inj_iw = [1,3,8]
            prod_iw = [2,4,5,6,7,9]
        elseif script == 2
            @info "Пятиточка, 4 нагнетательные ⁘ обр."
            inj_iw = [2,4,6,8]
            prod_iw = [1,3,5,7,9]
        elseif script == 3
            @info "Пятиточка, 5 нагнетательные ⁙ обр."
            inj_iw = [1,3,5,7,9]
            prod_iw = [2,4,6,8]
            #⁜
        elseif script == 4
            @info "Смещённая девятиточка +"
            inj_iw = [1,3,7,9]
            prod_iw = [2,4,5,6,8]
            #⁜
        else
            @info "Классическая девятиточка"
            inj_iw = [5]
            prod_iw = [1,2,3,4,6,7,8,9]
        end
    elseif nw == 36
        @info "36 скважин (модель от Айдара)"
        if script == 4
            @info "4 девятиточки от Айдара"
            inj_iw = [1,10,19,28]
            prod_iw = setdiff(1:36, inj_iw)
        end
    end

    return prod_iw, inj_iw
end
